#include "constant-advection.hh"
#include <math.h>
namespace Dune {
  namespace Initial {
    ConstantAdvection::ConstantAdvection(Dune::ParameterTree& ptree){
      set_initcond_para("dim",1.0);
      set_initcond_para("u0",0.0);
      set_initcond_para("v0",0.0);
      set_initcond_para("w0",0.0);

      //
      set_initcond_para("rectangle",0.0);
      set_initcond_para("circle",0.0);
      //
      set_initcond_para("wid1",0.0);
      set_initcond_para("wid2",0.0);
      set_initcond_para("wid3",0.0);
      set_initcond_para("r",0.0);
      set_initcond_para("cen1",0.0);
      set_initcond_para("cen2",0.0);
      set_initcond_para("cen3",0.0);
      //
      set_initcond_paras(ptree);
    }

    ConstantAdvection::~ConstantAdvection(){
    }

    void ConstantAdvection::set_initcond_paras(Dune::ParameterTree& ptree){
      para["dim"]=ptree.get("grid.dimension",0.0);
      para["u0"]=ptree.get("equation.u0",0.0);
      para["v0"]=ptree.get("equation.v0",0.0);
      para["w0"]=ptree.get("equation.w0",0.0);
      para["rectangle"] = ptree.get("initial.rectangle",0.0);
      para["circle"] = ptree.get("initial.circle",0.0);
      para["cen1"]=ptree.get("initial.cen1",0.0);
      para["cen2"]=ptree.get("initial.cen2",0.0);
      para["cen3"]=ptree.get("initial.cen3",0.0);
      para["r"]=ptree.get("initial.r",0.0);
      para["wid1"]=ptree.get("initial.wid1",0.0);
      para["wid2"]=ptree.get("initial.wid2",0.0);
      para["wid3"]=ptree.get("initial.wid3",0.0);
    }

    void ConstantAdvection::initial_condition(double* location, double* holder){
      if (para["rectangle"]==1.0)
        rectangle(location,holder);
      if (para["circle"]==1.0)
        circle(location,holder);
      return;
    }

    void ConstantAdvection::rectangle(double* location, double* holder){
      int dim = (int)para["dim"];
      int nvar = (int)para["nvar"];
      double cen[dim], wid[dim];
      if(dim>=1){
        cen[0]=para["cen1"];
        wid[0]=para["wid1"];
      }
      if(dim>=2){
        cen[1]=para["cen2"];
        wid[1]=para["wid2"];
      }
      if(dim==3){
        cen[2]=para["cen3"];
        wid[2]=para["wid3"];
      }
      // holder
      holder[0]=1.0;
      for(int i=0; i<dim; i++)
        holder[0]*=(std::abs(location[i]-cen[i])<wid[i]? 1.0 : 0.0);
      for(int i=1; i<nvar; i++)
        holder[i]=0.5;
      return;
    }

    void ConstantAdvection::circle(double* location, double* holder){
      int dim = (int)para["dim"];
      int nvar = (int)para["nvar"];
      double cen[dim], radius;
      if(dim>=1)
        cen[0]=para["cen1"];
      if(dim>=2)
        cen[1]=para["cen2"];
      if(dim==3)
        cen[2]=para["cen3"];
      radius=para["r"];
      // distance
      double distance=0.0;
      for(int i=0; i<dim; i++)
        distance=distance+std::pow(location[i]-cen[i],2);
      distance=std::sqrt(distance);
      // holder
      holder[0]=(distance<radius? 1.0 : 0.0);
      for(int i=1; i<nvar; i++)
        holder[i]=0.5;
      return;
    }
 } // Initial
}  // Dune
